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ABSTRACT 


Increasingly,  nuclear  explosion  monitoring  is  focusing  on  detection,  location,  and  identification  of  small 
events  recorded  at  regional  distances.  Because  Earth  structure  is  highly  variable  on  regional  scales,  locating 
events  accurately  at  these  distances  requires  the  use  of  region-specific  models  to  provide  accurate  travel 
times.  Improved  results  have  been  achieved  with  composites  of  ID  models  and  with  approximate  3D 
models  with  simplified  upper  mantle  structures,  but  both  approaches  introduce  nonphysical  boundaries  that 
are  problematic  for  operational  monitoring  use.  Ultimately,  what  is  needed  is  a  true,  seamless  3D  model  of 
the  Earth. 

Towards  that  goal,  we  have  developed  a  3D  tomographic  model  of  the  P  velocity  of  the  crust  and  upper 
mantle  for  the  region  of  southcentral  Asia  centered  around  the  Tibetan  Plateau.  Our  model  is  derived  from 
almost  140,000  Pn  picks  for  more  than  5300  events  recorded  at  563  stations  from  a  Ground  Truth  (GT) 
dataset  assembled  by  Los  Alamos  National  Laboratory  (LANL).  Our  starting  model  is  the  a  priori  model  of 
East  Asia  developed  by  LANL,  which  is  based  on  various  global  and  regional  studies.  The  topmost  layers 
come  from  the  Laske  and  Masters  global  sedimentary  model  from  1997.  As  our  dataset  lacks  the  resolution 
to  improve  this  sedimentary  portion  of  our  model,  we  fix  the  velocity  and  depth  of  these  layers,  as  well  as 
the  depths  of  the  major  mantle  discontinuities  (Moho,  410  km,  660  km).  We  invert  for  P  velocities  from  the 
crust  down  through  the  upper  mantle,  along  with  source  and  receiver  terms  to  account  for  the  effects  of 
event  mislocation  and  fine-scale  structure  near  the  receiver  not  accounted  for  in  the  crustal  model.  Eorward 
calculations  are  made  using  our  own  implementation  of  the  Um  and  Thurber  ray  pseudo-bending  approach 
(Ballard,  2008,  these  Proceedings)  with  full  enforcement  of  Snell’s  Law  in  3D  at  the  major  discontinuities. 
Because  large  portions  of  the  model  are  underconstrained,  we  apply  both  damping  and  regularization  to 
force  the  inversion  to  update  the  initial  model  only  where  good  data  coverage  is  available.  We  validate  the 
model  by  performing  several  inversions  with  random  portions  of  the  dataset  omitted  and  then  testing  the 
predictive  capability  of  the  model  against  those  portions  compared  with  AK135.  We  test  the  location 
performance  of  the  model  by  relocating  the  GT  events  using  our  model  and  using  AK135. 
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OBJECTIVES 

Our  goal  is  to  better  predict  Pn  travel  times  for  an  area  of  south  central  Asia  centered  around  the  Tibetan 
Plateau.  This  is  an  area  of  high  monitoring  interest,  as  it  encompasses  the  Indian,  Pakistani,  and  Chinese 
test  sites,  as  well  as  the  former  Soviet  test  site  at  Shagan  River  (now  in  Kazakhstan).  It  also  spans  several 
other  countries  that  could  be  candidates  for  nuclear  testing.  Thus,  locating  events  accurately  in  this  region 
is  tremendously  important  for  tracking  treaty  compliance  and  nuclear  weapon  proliferation.  To  accomplish 
our  goal,  we  have  developed  a  three-dimensional  (3D)  P  velocity  model  of  the  crust  and  upper  mantle  in 
this  region.  This  is  one  of  the  most  tectonically  complex  areas  on  the  Earth,  due  to  the  collision  of  the 
Indian  and  Asian  continents,  which  has  led  to  the  highest  mountain  ranges  in  the  world,  as  well  as  the 
thickest  crust  and  lithosphere. 

RESEARCH  ACCOMPLISHED 

Dataset 

Our  data  are  drawn  from  an  East  Asia  GT  dataset  collected  over  the  past  several  years  by  Ground-Based 
Nuclear  Explosion  Monitoring  Research  and  Development  (GNEMRD)  researchers  at  LANE.  Considerable 
effort  has  gone  into  the  collection  of  this  dataset  to  draw  on  all  possible  sources  of  event  catalogs  for  this 
region,  and  then  to  carefully  screen  them.  Events  are  drawn  from  global  catalogs  from  the  NEIC  (National 
Earthquake  Information  Center),  and  ISC  (International  Seismological  Centre),  as  well  as  from  various 
regional  catalogs  including  the  Annual  Bulletin  of  Chinese  Earthquakes  (ABCE)  (Lee  et  al.,  2002).  All 
events  in  the  LANE  GT  set  are  assigned  a  GT  level  based  on  the  Bondar  criteria  (Bondar  et  al.,  2004),  and 
for  this  project,  only  events  with  GT  levels  of  25  km  or  better  were  used.  For  each  of  these  events,  picks 
from  the  various  catalogs  were  merged,  removing  redundant  and/or  outlier  arrivals  based  on  a  preferred 
catalog  hierarchy.  Events  without  a  priori  GT  information  were  relocated  with  the  merged  arrival  set, 
comparing  the  results  to  the  Bondar  criteria  (Begnaud  et  al.,  2006).  To  focus  on  P  velocity  in  the  upper 
mantle,  we  chose  to  use  only  Pn  phases  and  have  limited  source  to  receiver  path  lengths  to  be  25  degrees  or 
less.  This  resulted  in  approximately  1 10,000  Pn  picks  for  more  than  5300  events  recorded  at  563  stations 
(Figure  1). 


5389  events,  3156  stations 


Figure  1.  Map  showing  the  seismic  events  (red  circles)  and  stations  (green  triangles)  used  in  this 

study.  The  magenta  dashed  line  shows  the  extent  of  the  tomographic  model  we  invert  for, 
while  the  lines  of  hlack  circles  show  the  positions  of  cross  sections  through  our  model 
presented  later. 
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Model  Structure 

One  of  the  areas  in  which  we  believe  our  model  is  an  improvement  over  many  of  the  previous  tomographic 
models  of  the  region  is  in  the  structure  of  the  model.  Many  previous  studies  have  used  geographic  grids  that 
are  regularly  spaced  in  latitude  and  longitude.  These  types  of  grids  are  easy  to  construct  and  to  use  for 
interpolation  (locating  neighbor  node  positions  within  them  is  trivial),  but  they  are  inefficient  because  they 
compress  northward  as  the  longitude  meridians  converge.  We  instead  use  a  grid  that  is  truly  equally  spaced, 
regardless  of  latitude,  as  is  shown  in  Figure  2.  This  grid  is  derived  by  progressive  subdivision  of  an 
icosahedron  with  12  nodes  to  a  set  of  nodes  with  the  desired  spacing  (see  Ballard  et  al.,  2008a,  these 
Proceedings,  for  more  details).  As  with  a  regular  geographic  grid,  values  at  intermediate  positions  are 
found  by  interpolation  of  values  at  the  neighbor  nodes,  but  for  this  type  of  grid  the  neighbor  nodes  are 
found  by  a  standard  walking  triangle  search.  To  avoid  the  complications  of  dealing  with  ray  paths  that  cross 
in  and  out  of  our  model,  we  use  a  global  grid,  but  the  grid  spacing  is  variable,  being  1  degree  in  Eurasia  and 
2  degrees  everywhere  else  in  the  world.  The  radial  structure  at  each  of  these  geographic  node  positions  is 
represented  with  a  standard  set  of  layers  that  extend  from  the  surface  of  the  Earth  to  the  center  of  the  Earth: 
sediments,  upper  crust,  mid-crust,  lower  crust,  uppermost  mantle,  etc.  The  thicknesses  (layers  can  pinch 
out),  top  and  bottom  depths,  and  top  and  bottom  velocities  of  each  of  these  layers  can  all  change  with 
geographic  position,  though  we  do  assume  that  the  velocity  within  a  layer  is  fit  by  a  linear  gradient  between 
the  top  and  bottom  values.  Values  at  intermediate  positions  within  layers  are  interpolated.  Our  model 
structure  is  described  in  more  detail  in  the  companion  presentation  by  Ballard  et  al.  (2008a,  these 
Proceedings). 


Figure  2.  Geographic  grid  for  the  tomographic  model.  Color  is  proportional  to  topography  (red  is 

high,  green  is  medium,  and  blue  is  low).  Each  node  marks  the  position  of  a  layered  profile 
extending  from  the  surface  to  the  center  of  the  Earth.  The  hlack  line  around  the  Indian 
subcontinent  and  the  Tibetan  Plateau  shows  the  portion  of  the  overall  model  that  will  be 
updated  by  our  tomographic  inversion. 
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Starting  Model 

As  mentioned  above,  even  though  our  goal  is  improved  structure  for  a  limited  region  of  southern  Asia,  we 
have  found  it  to  be  more  straightforward  to  work  with  a  global  model.  For  a  background,  we  use  the  crustal 
structure  from  the  Laske,  Master,  and  Reif  Crust  2.0  model  (Laske  et  al.,  2008)  on  top  of  an  AK135  mantle. 
For  Europe,  North  Africa,  and  Asia,  we  use  the  a  priori  GNEMRD  Unified  Model  created  collaboratively 
by  LANE  and  Lawrence  Livermore  National  Laboratory  (LLNL)  (Begnaud  et  al.,  2004;  Pasyanos  et  al., 
2004).  This  model  has  been  shown  to  provide  better  Pn  travel  times  than  AK135  for  Eurasia  and  North 
Africa  (Elanagan  et  al.,  2007),  but  not  for  southern  Asia  (Young  et  al.,  2008),  hence  our  motivation  for  this 
study.  Shallow  sedimentary  layers  are  taken  from  Crust  2.0,  while  the  deeper  crustal  layers  are  taken  from 
various  focused  regional  studies.  The  model  also  includes  upper  mantle  structure  from  the  3SMAC  model 
(Nataf  and  Ricard,  2004),  but  we  choose  not  to  use  this  because  it  is  not  based  on  seismic  data.  Instead,  we 
use  AK135  mantle  structure. 

Limitations  on  the  Inversion 

Based  on  the  ray  paths  connecting  our  sources  and  receivers,  we  constructed  the  polygon  shown  in 
Figures  1  and  2  to  limit  the  extent  of  our  inversion.  Only  the  velocities  at  the  nodes  within  this  polygon 
were  updated.  Further,  we  do  not  allow  either  the  thicknesses  or  the  depths  of  layers  to  change:  only  P 
velocities  at  the  tops  and  bottoms  of  layers  can  be  modified.  We  only  include  ray  paths  in  the  inversion  that 
have  3/4  or  more  of  their  path  lengths  lying  within  our  polygon.  For  such  rays,  observed  times  are  always 
reduced  by  the  travel  time  for  the  portion  of  the  path  external  to  the  polygon  before  inversion.  Because  we 
are  using  Pn  data  only,  we  limit  our  inversion  to  rays  that  turn  below  the  Moho  but  above  the  410  km 
discontinuity.  Hence,  our  inversion  is  limited  to  the  crust  and  the  portion  of  the  upper  mantle  above  the 
transition  zone. 

Inversion  Formulation 

Our  inversion  minimizes  the  travel  time  residuals  formed  by  subtracting  predicted  Pn  travel  times  from 
observed  travel  times.  We  assume  that  the  travel  time  residual  along  a  path  from  a  source  i  to  a  receiver) 
can  be  calculated  by  multiplying  the  required  slowness  perturbations  along  the  path  by  the  lengths  of  the 
segments  along  the  path,  i.e., 

where  d,,  is  the  length  of  a  path  along  the  1th  segment  and  is  the  slowness  perturbation  along  that 

segment.  With  perfect  data  and  perfect  coverage  of  the  entire  model,  this  equation  could  be  inverted 
directly  to  solve  for  the  slowness  perturbations  (and  hence  the  velocity  structure).  However,  the  data  are 
never  perfect,  nor  is  the  coverage,  so  a  best  model  is  solved  for  iteratively  using  a  least-squares 
minimization  of  the  residuals. 

In  tomographic  inversions  where  the  model  is  structured  as  a  set  of  blocks,  it  is  the  path  length  through 
those  blocks  that  provide  the  weights  to  be  multiplied  by  the  slowness  perturbations  to  get  the  residuals. 
Because  our  model  is  represented  by  a  set  of  nodes  instead  of  blocks,  the  formulation  is  slightly  different. 
We  divide  each  path  length  into  equal  length  subsegments  10  km  long,  then  for  each  subsegment,  we 
interpolate  from  the  surrounding  nodes  to  get  the  appropriate  slowness  to  multiply  times  the  subsegment 
length.  Hence,  in  our  formulation,  we  are  solving  for  the  slowness  perturbation  at  the  nodes,  and  the 
coefficients  multiplying  each  of  the  node  slowness  perturbations  are  a  combination  of  the  path  subsegment 
lengths  and  the  weights  of  the  nodes  used  for  the  interpolation.  Because  our  data  have  noise,  we  cannot  find 
a  model  that  gives  a  zero  residual  everywhere  but  instead  seek  to  minimize  the  squared  residual: 

£  =||AA5' -  A/||^ ,  (2) 

where  A  is  an  «  (number  of  observations)  by  m  (number  of  nodes)  matrix  of  node  coefficients.  As  is  an  m 
by  1  array  of  slowness  perturbations,  and  At  is  an  n  by  1  array  of  travel-time  residuals. 
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Because  portions  of  the  model  are  under-constrained,  we  impose  damping.  In  addition,  we  add  horizontal 
and  vertical  smoothing  (regularization)  to  produce  models  that  are  consistent  with  our  geophysical 
intuition.  Thus,  the  overall  equation  we  invert  is  of  this  form: 

s  =  II AAv  -  Af||^  +  A||A,s'||^  +  y  ||LAj'||^  ,  (3) 

where  A  is  the  damping  coefficient,  y  is  the  smoothing  factor,  and  L  is  the  matrix  of  smoothing  operators. 
To  control  the  slowness  adjustments  in  regions  of  the  model  where  sampling  is  very  poor,  we  use  separate 
damping  coefficients  for  each  node  slowness  perturbation  and  make  the  damping  coefficient  for  each 
inversely  proportional  to  the  ray  coefficients,  i.e.,  nodes  with  very  little  path  sampling  are  highly  damped, 
while  those  with  high  path  sampling  are  damped  very  little. 

Travel- Time  Prediction  Calculation  and  Ray  Paths 

Another  area  in  which  we  believe  this  study  improves  over  previous  tomographic  studies  of  the  area  is  in 
our  calculation  of  travel  times  and  ray  paths.  For  each  iterative  model,  beginning  with  the  starting  model, 
we  calculate  true  3D  ray  paths  using  our  own  implementation  of  the  Um  and  Thurber  (1987) 
pseudo-bending  methodology.  Our  implementation  requires  rays  to  honor  Snell’s  Law  at  discontinuities 
(in  this  case  the  Moho),  and  it  also  avoids  local  minima  by  solving  for  a  ray  (refracted,  diffracted,  or 
reflected)  in  each  layer,  then  finding  the  first  refracting  arrival.  We  have  shown  our  implementation  to  be 
highly  accurate  (Young  et  al.,  2008),  but  also  computationally  expensive,  so  we  have  designed  a  distributed 
computing  framework  to  allow  a  set  of  rays  to  be  run  simultaneously  on  multiple  processors.  For  our 
calculations  for  this  model,  we  used  a  set  of  nine  16-processor  PC-architecture  machines  for  a  total  of  144 
processors.  For  more  on  the  ray  bending  algorithm  implementation,  we  refer  the  reader  to  the  companion 
paper  by  Ballard  et  al.  (2008b). 

Travel- Time  Prediction  Improvement 


The  ultimate  goal  of  this  project  is  to  develop  a  model  that  can  provide  better  Pn  travel  time  predictions  for 
our  area.  Our  model  shows  clear  improvements  over  the  AK135  model  for  our  dataset.  Figure  3  shows 
mean  and  standard  deviations  of  the  observed  minus  prediction  residuals  for  both  models. 


Distance  (degrees) 


Figure  3.  (left)  One  degree  station  to  event  distance  bins  of  mean  plus  and  minus  one  standard 
deviation  for  starting  model  (red)  and  final  model,  (right)  One  degree  station  to  event 
distance  bins  of  standard  deviation. 
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The  means  of  the  residuals  at  all  distances  are  negative  for  AK135,  indicating  that  the  data  require  a  faster 
model.  Further,  the  standard  deviations  are  large  (over  2  seconds),  suggesting  that  there  is  considerable 
lateral  heterogeneity  in  the  region. 

For  our  model,  the  overall  RMS  of  the  residuals  has  been  reduced  from  3.2  seconds  to  2.5  seconds,  an 
improvement  of  22%.  This  suggests  that  we  are  modeling  a  substantial  portion  of  the  lateral  heterogeneity 
in  the  region.  The  bulk  of  the  improvement  comes  beyond  8  degrees,  suggesting  that  it  is  the  changes  in  the 
upper  mantle  portion  of  our  model  that  are  the  most  important.  Our  model  decreases  the  bias,  especially 
beyond  8  degrees,  but  does  not  remove  it.  We  believe  that  the  reason  for  this  is  that  our  damping 
coefficients  were  too  strong,  hence  the  model  was  unable  to  move  far  enough  away  from  the  starting 
model.  In  subsequent  versions  of  the  model,  we  plan  to  vary  the  damping  coefficients  to  see  if  we  can 
remove  the  bias  entirely. 

Discussion  of  the  Tomographic  Model 

Figure  4  shows  percent  velocity  changes  from  the  staring  model  at  the  bottom  of  the  lowermost  crustal 
layer,  at  the  top  of  the  uppermost  mantle  layer,  at  100  km,  and  at  200  km.  To  help  interpret  these  pictures, 
at  each  depth  we  also  show  a  map  of  the  node  ray  hit  count  at  each  depth.  Generally  ray  hit  counts  are  good 
across  the  model  for  all  depths,  but  there  are  lower  hit  counts  at  the  edges  of  the  model,  especially  in  the 
crust  where  the  ray  paths  are  much  more  vertical  than  in  the  mantle. 

Beginning  with  the  shallowest  structure  and  working  downwards,  we  can  see  that  overall  the  adjustments  to 
the  velocity  at  the  base  of  the  crust  are  much  less  than  in  the  mantle.  This  suggests  that  the  crustal  portion 
of  our  starting  model  is  doing  a  good  job  of  fitting  the  data,  as  expected.  Still,  our  results  suggest  that  the 
crustal  velocities  beneath  the  Tibetan  Plateau  are  lower  than  in  the  starting  model,  and  that  the  velocities  to 
the  NW  and  SE  of  the  plateau  need  to  be  higher  to  match  the  observations. 

The  strongest  velocity  changes  are  seen  at  100-km  depth,  and  they  are  predominantly  velocity  increases, 
which  is  consistent  with  the  negative  bias  of  the  starting  model  seen  in  Figure  3.  Velocities  have  been 
increased  by  2%  and  greater  throughout  most  of  the  region,  with  the  exception  of  low-velocity  regions  near 
the  Tarim  Basin  and  in  the  southeast  part  of  China.  Both  the  positive  and  negative  features  seem  to  extend 
upward  to  the  Moho  and  downward  to  200  km,  though  not  as  strongly.  Depth  extent  can  be  better  assessed 
with  the  vertical  cross  sections  shown  in  Figure  5  for  the  profiles  shown  in  Figure  1.  We  can  see  a 
high-velocity  region  extending  from  the  Indian  subcontinent  underneath  the  Tibetan  Plateau,  then  being 
sharply  truncated  by  low-velocity  material  at  the  Tarim  Basin  for  slices  that  cross  that  basin  farther  to  the 
east.  Both  these  high  and  low-velocity  anomalies  extend  from  the  top  of  the  mantle  down  to  at  least 
200  km.  How  much  deeper  they  extend  is  not  clear  because  of  our  poor  ray  coverage  at  depth. 
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Figure  4.  (left  side,  top  to  bottom)  LoglO  ray  hit  counts  for  nodes  at  the  base  of  the  crust,  the  top  of 
the  mantle,  the  top  of  the  100-km-depth  layer,  and  the  top  of  the  200-km-depth  layer, 
(right  side,  top  to  bottom)  Percent  P  velocity  change  from  starting  model  for  the  same 
nodes  and  depths. 
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%  Vp  change,  slice  #1  Lat/Lon  16.0/65.3  to  47.8/80.2 
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Figure  5.  (left  side,  top  to  bottom)  LoglO  ray  hit  counts  for  cross  section  slices  shown  in  Figure  1. 

(right  side,  top  to  bottom)  Percent  P  velocity  change  from  the  starting  model  for  the  same 
cross  sections. 


CONCLUSIONS  AND  RECOMMENDATIONS 


We  have  developed  a  3D  model  of  the  P  velocity  structure  of  the  crust  and  upper  mantle  for  a  limited 
region  of  south  central  Asia.  Our  approach  differs  from  previous  efforts  in  several  important  ways.  First, 
we  utilize  only  high-quality  GT  data,  collected  hy  LANE  GNEMRD  researchers  over  the  past  several 
years.  Second,  our  model  is  represented  geographically  with  a  grid  of  truly  equally  spaced  nodes,  and 
radially  with  a  set  of  standard  layers  with  variable  properties.  Third,  we  trace  rays  and  calculate  travel  times 
using  a  3D  pseudo-hending  approach  that  enforces  Snell’s  law  across  discontinuities.  The  resulting  model 
provides  a  22%  improvement  in  RMS  travel-time  residual  over  our  starting  model,  which  is  taken  from  the 
LANL/LLNL  a  priori  unified  geophysical  model  that  encompasses  our  region.  Overall,  our  model  is  faster 
than  the  starting  model,  hence  reducing  the  negative  bias  in  residuals,  but  our  model  still  has  a  negative 
bias,  probably  because  damping  is  too  strong.  In  future  versions  of  the  model,  we  plan  to  experiment  with 
changing  damping  and  regularization  parameters  to  improve  our  model  further. 
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